function P=plya(A,B,g)
%Solve the unique positive definite solution to the Parametric Riccati matrix equation 
%
%                               A'P+PA-PBB'P=-gP
%
%via symbolic computation.
%Please input A and B defined by sym function.

%Bin Zhou, Harbin Institute of Technology, binzhou@hit.edu.cn
[n,m]=size(B);
C=kron(A+g/2*eye(n),eye(n))+kron(eye(n),A+g/2*eye(n));
B1=B*B.';
b=B1(:);
% b=vec(B*B.');
W=simplify(inv(C)*b);
W=invvec(W);
P=inv(W);
P=simplify(P);